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Abstract 

We present the first results from BlackHat, an automated C++ program for calculating one-loop 
qq \ amplitudes. The program implements the unitarity method and on-shell recursion to construct 



amplitudes. As input to the calculation, it uses compact analytic formulas for tree amplitudes 
for four-dimensional helicity states. The program performs all related computations numerically. 



CO 
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We make use of recently developed on-shell methods for evaluating coefficients of loop integrals, 



introducing a discrete Fourier projection as a means of improving efficiency and numerical stability. 
We illustrate the numerical stability of our approach by computing and analyzing six-, seven- and 
eight-gluon amplitudes in QCD and comparing against previously-obtained analytic results. 

PACS numbers: ll.15.Bt, 11.55.Bq, 12.38.Bx 
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I. INTRODUCTION 



The Large Hadron Collider (LHC) will soon begin exploration of the electroweak sym- 
metry breaking scale. It is widely anticipated that physics beyond the Standard Model will 
emerge at this scale, leading to a breakthrough in our understanding of TeV-scale physics. 
A key ingredient in this quest is the precise understanding of the expected Standard Model 
backgrounds to new physics from both electroweak and QCD processes. In the absence of 
such an understanding, new physics signals may remain hidden, or backgrounds may be 
falsely identified as exciting new physics signals. 

Quantitatively reliable QCD predictions require next-to-leading order (NLO) calcula- 
tions [1]. For a few benchmark processes, such as the rapidity distribution of electroweak 
vector bosons 3], the transverse-momentum distribution of the Z boson at moderate pr> 
and the total cross sections for production of top quark pairs and of Higgs bosons J3( , the 
higher precision of next-to-next-to- leading order (NNLO) results may be required 1|. For 
most other processes, NLO precision should suffice. However, there are many relevant pro- 
cesses that need to be computed, particularly those with high final-state multiplicity. Such 
processes are backgrounds to the production of new particles that have multi-body decays. 
To date, no complete NLO QCD calculation involving four or more final-state objects (par- 
ticles or jets) is available. (In electroweak theory, however, e + e~ — > 4 fermions has been 
evaluated |4| using the integral reduction scheme of Denner and Dittmaier |5j.) NLO correc- 
tions require as ingredients both real-radiative corrections and virtual corrections to basic 
amplitudes. The structure of the real-radiative corrections — isolation of infrared singu- 
larities and their systematic cancellation against virtual-correction singularities — is well 
understood, and there are general methods for organizing them Indeed, the most 

popular of these methods, the Catani-Seymour dipole subtraction method [8j , has now been 
implemented in an automatic fashion [9]. The infrared divergences of virtual corrections, 
needed to cancel the divergences from integrating real radiation over phase space, are also 
understood in general jf], [lO]. The main bottleneck to NLO computations of processes with 
four or more final-state objects has been the evaluation of the remaining ingredients, the 
infrared-finite parts of the one-loop virtual corrections. 

As the number of external particles increases, the computational difficulty of loop- 
amplitude calculations using traditional Feynman diagrams grows rapidly. Technologies 
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that have proven useful at tree level, such as the spinor-helicity formalism 11], do not suf- 
fice to tame these difficulties. In the past few years, several classes of new methods have 
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been proposed to cope wit 
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which are based on t 
tude must satisfy [40 



le analytic properties of unitarity and factorization that any ampli- 



411 ]. These methods are efficient, and display very mild growth in 



required computer time with increasing number of external particles, compared to a tradi- 
tional Feynman-diagrammatic approach. The improved efficiency emerges from effectively 
reducing loop calculations to tree-like calculations. Efficient algorithms can then be em- 
ployed for the tree-amplitude ingredients. 

oped in 
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42] 



One of the principal on-shell technologies is the unitarity method, originally deve 
calculations of supersymmetric amplitudes with more than four 1 external particles 
An early version combining unitarity with factorization properties was used to compute 
the one-loop amplitudes for e + e~ — > Z — > 4 partons and (by crossing) for amplitudes 
entering pp — > W, Z + 2 jets 181 ] . (The latter have been incorporated into the NLO program 
MCFM |43|] . ) This calculation introduced the concept of generalized unitarity 4l|] as an efficient 
means for performing loop computations. It improves upon basic unitarity because it isolates 
small sets of terms, and hence makes use of simpler on-shell amplitudes as basic building 
blocks. On-shell methods have already led to a host of new results at one loop, including the 



computation of non-trivial amplitudes in QCD with an arbitrary number of external legs 25, 



26 



27 



28 



44j. This computation goes well beyond the scope of traditional diagrammatic 



computations, and provides a clear demonstration of the pow er of the methods. The reader 
may find recent reviews and further references in refs. jll. \s<\. 

The next challenge is to move beyond analytic calculations of specific processes or 
classes of processes to produce a complete, numerically stable, efficient computer code 
based on these new developments. Here we report on an automated computer program 
- BlackHat - - based on on-shell methods, with the stability and efficiency required to 
compute experimentally-relevant cross sections. Other researchers are constructing numeri- 



cal programs 



35, 



36 



37 



38. 



391 ] based on related methods 



31 



37j, 



39|. 



1 The earlier dispersion relation approach 40( had not been used to construct amplitudes with more than 
two kinematic invariants. 



3 



On-shell methods rely on the unitarity of the theory [40j and on its factorization proper- 
ties, which together require that the poles and branch cuts of amplitudes correspond to the 
physical propagation of particles. In general, any one-loop amplitude computed in a quantum 
field theory contains terms with branch cuts, and also purely rational terms, that is, terms 
that have no branch cuts and are rational functions of the external momentum invariants (or 
more precisely of spinor products). The cut- containing pieces can be determined from uni- 



tarity cuts, in which the intermediate states may be treated four-dimensionally [17|,|42|]. Only 
products of tree-level, four- dimensional helicity amplitudes are needed for this step. The ra- 
tional terms have their origin in the difference between D = 4 — 2e and four dimensions when 
using dimensional regularization. They can be obtained 2 within the unitarity metho d by 
keeping the full D-dimensional dependence of the tree amplitudes 
Alternatively, to obtain the rational terms, one can use on-shell recursion 
struct the rational remainder from the loop amplitudes' factorization poles [26j, |28|, |44( . We 
will follow the latter route in this paper. 

A generic one-loop amplitude can be expressed in terms of a set of scalar master in- 
tegrals m ultiplied by various rational coefficients, along with the additional purely rational 

The relevant master integrals depend on the masses of the physical 




terms 



46, 



43, 



18 



states that appear, but otherwise require no process-specific computation. At one loop, they 



consist of box, triangle, bubble anc 



integrals are known analytically 51 



for massive particles) tadpole integrals. The required 



52] 



Our task is therefore to determine the coefficients in front of these integrals for each 
process and helicity configuration. We do so using generalized cuts 

HflHQ- Britto, 

Cachazo and Feng (BCF) observed 2jj that with complex momenta one can use quadruple 
cuts to solve for all box coefficients, because massless three-point amplitudes isolated by cuts 
do not vanish as they would for real massless momenta. Moreover, the solution is purely 
algebraic, because the loop momentum of the four-dimensional integral is completely frozen 
by the four cut conditions, and a given quadruple cut isolates a unique box coefficient. This 
provides an extremely simple method for computing box-integral coefficients. Continuing 
along these lines, Britto, Buchbinder, Cachazo, Feng, and Mastrolia have developed efficient 



2 This fact is closely connected to van Neerven's important observation that dispersion relations for Feynman 
integrals converge in dimensional regularization 45 1. 
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analytic techniques [29j for evaluating generic one-loop unitarity cuts to compute triangle 
and bubble coefficients. They use spinor variables and compute integral coefficients via 
residue extraction. 

For the purposes of constructing a numerical code, we use a somewhat different ap- 
proach. For triangle integrals, we can impose at most three cut conditions. This leaves a 
one-parameter family of solutions. These conditions no longer isolate the triangle integral 
uniquely, as a number of box integrals will share the same triple cut. Similar considerations 
apply to the ordinary two-particle cuts needed to obtain bubble coefficients. As discussed 
by Ossola, Papadopoulos and Pittau (OPP) [311 ] . one can construct a general parametric 
form for the integrand. This form can be understood as a decomposition of the loop mo- 
mentum in terms of components in the hyperplane of external momenta and components 
perpendicular to this hyperplane [35)]. Coefficients of the various master integrals can be 



extracted by comparing the expressions obtained from Feynman graphs with the general 
parametric form, using values of the loop momentum in which different combinations of 
propagators go on shell. For the quadruple cut, this leads to a computation identical to the 
method of ref. 21] once one further replaces sums of Feynman diagrams by tree amplitudes. 
OPP solve the problems of box contributions to triangle coefficients, and of box and triangle 
contributions to bubble coefficients, iteratively by subtracting off previously-determined con- 
tributions and solving a particular system of equations numerically. In the OPP approach, 
the rational terms can be determined by keeping the full D-dimensional dependence in all 



terms 



31 



38]. 



Forde's alternative approach makes use of a complex-valued parametrization of the loop 



momenta [34|J (similar to the one used in refs. [13|, |31|) and exploits the different func- 
tional dependence on the complex parameters to separate integral contributions to a given 
triple or ordinary cut. We develop this method one step further, and introduce a discrete 
Fourier projection in these complex parameters, in conjunction with an OPP subtraction of 
previously-determined contributions 3JJ. The projection isolates the desired integral coef- 
ficients efficiently, while maintaining good numerical stability in all regions of phase space. 



35 



38 



391 ] from solving a system of linear 



It minimizes the instabilities that may arise 
equations in regions where the system degenerates. 

We compute the rational remainder terms using loop-level on-she 



tions 



26 



28 



44j, analogous to the recursion relations at tree level 23 



1 recursion rela- 



24] developed by 
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Britto, Cachazo, Feng and Witten (BCFW). At tree level, gauge-theory amplitudes can be 
constructed recursively from lower-point amplitudes, by applying a complex deformation to 
the momenta of a pair of external legs, keeping both legs on shell and preserving momentum 
conservation. The proof relies only on the factorization properties of the theory and on 
Cauchy's theorem, so the method can be applied to a wide variety of theories. At loop level, 
the construction of an analogous recursion relation for the rational terms requires addressing 
a number of subtleties, including the presence of spurious singularities. These issues can 
and have been addressed for specific infinite series of one-loop helicity amplitudes, allowing 



their recursive construction |26|, [27|, |28|, |44j . In order to more easily automate the method of 
ref. 26|, we modify how the spurious singularities are treated, making use of the availability 



of the integral coefficients within the numerical program, in a manner to be described below. 

In any numerical method, the finite precision of a computation means that instabilities can 
arise, occasionally leading to substantial errors in evaluating an amplitude at a given point 
in phase space. We introduce simple tests for the stability of the evaluation. Principally, we 
check that the sum of bubble integral coefficients agrees with its known value, and we check 
for the absence of spurious singularities in this sum. A comparison with known analytic 
answers for a variety of gluon amplitudes shows that these two tests suffice to detect almost 
all instabilities. If a test fails, we consider the point to be unstable. Various means of dealing 
with unstable points have been discussed jj a [ 



54 



55. 



561 ] . We simply re-evaluate the 

(. Doing 



fairly small fraction of unstable points at higher precision using the QD package [5 
so, we still have an average evaluation time of less than 120 ms for the most complicated of 
the six-gluon helicity amplitudes, and subtantially better times for the simpler ones. Higher- 
precision evaluation has also been used recently in ref. 36( to handle numerically unstable 
points. 

Although BlackHat is written in C++, for algorithm development and prototyping, we 
found it extremely useful to use symbolic languages such as Maple [58( and Mathematica [59] , 
and in particular the Mathematica implementation of the spinor-helicity formalism provided 
by the package S@M 60]. At present BlackHat computes multi-gluon loop amplitudes. Once 
we implement a wider class of processes in the same framework, we intend to release the 
code publicly. 

The present paper is organized as follows. In section [Til we discuss how we compute the 
coefficients of the various integral functions, and introduce the discrete Fourier projection. 
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FIG. 1: The basis of scalar integrals: (a) box, (b) triangle, (c) bubble, and (d) tadpole. Each corner 
can have one or more external momenta emerging from it. The tadpole integral (d) vanishes when 
all internal propagators are massless. 



In section UTTI we outline the calculation of the purely-rational terms, describing in particular 
our treatment of the spurious singularities. We also introduce our criteria for ensuring the 
numerical stability of the computed amplitude. We show results for a number of gluon 
amplitudes with up to eight external legs in section HV 
defer a number of technical details to a future paper 6JJ 



and summarize in section |Vj We 



II. INTEGRAL COEFFICIENTS FROM FOUR-DIMENSIONAL TREE AMPLI- 
TUDES 

We begin by dividing the dimensionally-regularized amplitude into cut-containing and 
rational parts. We evaluate the cut parts using the four-dimensional unitarity method 1?], 
331 ] . To extract the box- integral coefficients we use the observation of BCF that the quadruple 
cuts freeze the loop integration 21]. For triangle and bubble integrals we use key elements 
of both the OPP [31] and Forde {34] approaches. In addition, we introduce a discrete 
Fourier projection for extracting the integral coefficients. (Alternative on-shell methods for 



29 



obtaining the integral coefficients have been given in refs. 

As the first step, we separate an n-point amplitude A n into a cut part C n and a rational 
remainder R n , 

A n = C n + R n . (2.1) 



The cut part is given by a linear combination of scalar basis integrals (46 



47|, 



48 



49, 



50 



5l|, 



c n = j2 dJi + E c * J s + E h A + E a i p i 



(2.2) 
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FIG. 2: The quadruple cut used to determine the coefficients of the box integrals. The loop 
momenta, flowing clockwise, are constrained to satisfy on-shell conditions. The blobs at each 
corner represent tree amplitudes. The dashed lines indicate the cuts. The external momenta are 
all outgoing. 



The integrals are scalar box, triangle, bubble and tadpole integrals, illustrated 

in fig. [TJ For massless particles circulating in the loop, the tadpole integrals vanish in 
dimensional regularization. The integral coefficients dj, q, a,i are rational functions of 
spinor products and momentum invariants of the kinematic variables, and are independent 
of the dimensional regularization parameter e. The index i runs over all distinct integrals of 
each type. The rational terms R n are defined by setting all scalar integrals to zero, 3 



Rn — A 7 



If ^0 



(2.3) 



Alternatively, the rational terms can be absorbed into the integral coefficients by keeping 
their full dependence 4 on e. 

In this paper, we obtain the integral coefficients at e = by using the unitarity 



method with four-dimensional 
dimensional spinor techniques 



oop momenta. This method allows us to use powerful four- 



11 



631 ] to greatly simplify the tree amplitudes that serve as 



basic building blocks. We will instead obtain the rational terms R n from on-shell recur- 



sion 



23 



24 l26l . |27] , as explained in the subsequent section. 



3 All contributions from the scalar integrals in eq. (|2.2j) are part of C n , including all 1/e 2 and 1/e pole 
terms, tt 2 factors, and pieces arising from the order e° term in the scalar bubble integral. 

4 The e dependence leads only to rational contributions, because it arises from integrals with (— 2e) com- 
ponents of loop momenta in the numerator. Each such integral can be rewritten as the product of e with 
a higher-dimensional integral, which possesses at most a single, ultraviolet pole in e, whose residue must 
be rational 13. |62|. 



8 



A. Box Coefficients 



Consider first the coefficients of the box integrals. We obtain them from the quadruple cut 
shown in fig. [2J The cut propa gato rs correspond to the four propagators of the desired box 
coefficient. As observed in ref. 21(, if we take the loop momentum to be four-dimensional, 



then the four cut conditions, 

Z? = m?, % = 1,2,3,4, (2.4) 

match the number of components of the loop momentum, leading to a discrete sum over two 
solutions for U. The integration is effectively frozen. The m; are the masses of the particles 
in the cut propagators, which in this paper are taken to vanish. The coefficient of any box 
integral is then given in terms of a product of four tree amplitudes, 

di = ~E^> (2-5) 

Z <T=± 



ju _ a tree a tree a tree a tree 
a i — ^(l) A {2) A (3) ^(4) 



(2.6) 



where the sum runs over the two solutions to the on-shell conditions, labeled by "+" and 
"— ". The four tree amplitudes in eq. (12.61) correspond to the tree amplitudes at the four 
corners of the quadruple cut depicted in fig. [2j 



The generic solution for was found in ref. 



211 ]. Simpler forms can be found for 



particular, but still fairly general, kinematical cases. In this paper, we focus on the case of 
massless particles circulating in the loop. When in addition at least one external leg, say 
leg 1, of the box integral shown in fig. [2] is also massless, that is K\ = 0, the two solutions 
to the on-shell conditions (I2.4[) can be written in a remarkably simple form, 









1*) 


2(1 T 




2<1=F|# 2 # 4 |1 ± ) 



I's J " n/^l* Mix > ^ J - 2(1^|« 4 |1 ± ) • (2 - 7j 

As illustrated in fig. (2J the A"j are the external momenta of the corners of the box integral 
under consideration and and |l ± )are Weyl spinors corresponding to the massless mo- 
mentum Ki, in the notation of refs. [63J. In massless QCD this solution covers all helicity 
configurations for amplitudes with up to seven external quarks or gluons, and a large frac- 
tion of the box coefficients for more external partons. (This solution may also be found in 



Risager's Ph.D. thesis 



6j.) 




FIG. 3: (a) The triple cut and (b) the ordinary double cut used to determine the coefficients of the 
triangle and bubble integrals. The loop momenta flowing clockwise, are constrained to satisfy 
on-shell conditions. The external momenta are all outgoing. 

The solution (12.71) has the advantage of making it manifest that Gram determinants enter 
only as square roots, one for each power of loop momenta in the numerator of box integrals. 
Indeed, the Gram determinant is given by the product of the spinor-product strings in the 
denominators of eq. (12.71) . 

A 4 = -2 $ 2 $ A |1+) (1+| $ 2 $ A |l-> , (2.8) 

where A 4 = det(2JQ • Kj), i,j = 1,2,3, is the box Gram determinant for K\ = 0. This 
property reduces the severity of numerical round-off error due to cancellations between 
different terms, in the regions of phase space where the Gram determinant vanishes. 



B. Triangle Coefficients from Discrete Fourier Projection 

To evaluate the coefficients of the triangle and bubble integrals, we make use of elements 
from the approaches of both OPP [31] and Forde [34]. First consider the triangle integrals. 
To obtain the coefficients q in eq. (12.21) we use the triple cut depicted in fig. [3]^a). In 
contrast to the quadruple cut, the triple cut does not freeze the integral, but leaves one 
degree of freedom which we denote by t. Moreover, the triple cut also contains box integral 
contributions. This makes the extraction of the triangle coefficients somewhat more intricate 
than the box coefficients. 

For massless internal particles, the solution of the cut condition If — (i = 1, 2, 3) 
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IS 



13, 



31 



24 1 



t ,~ 



Il(t) = K? + K£+-(iq;\>f\K;) + -{Ks\>f\Kr), 



2t 



(2.9) 



and, using momentum conservation, h{t) = h(t) — ifi, l 3 {t) = li(t)+K 3 . Here t is a complex 
parameter corresponding to the one component of the loop momentum not fixed by the cut 
condition. Following ref. 34J we have, 

,7#£ + SsKlt 



iff = 7a 



7 iff + Siif£ 



Kit 



-"fa 



[2.10) 



9 C 1 ci ' O / 9 CI CI 

Y - 6i6 3 y - SiS 3 

with S\ = Kf, S 3 = iff, and iff and K 3 are both massless. (In comparison with ref. [34j, we 
have rescaled and relabeled these massless momenta, and here we take all external momenta 
to be outgoing.) The variables a, a' and 7 are defined as follows, 

SziSi-Y , Si(s 3 ~y 



a 



SiS 3 - Y 



a: 



SiS 3 - Y 



7 



7± = -if 1 -if 3 ±v / A, (2.11) 



where 



A = - det(ifi • Kj) = {Ki ■ if 3 



K\K, 



3 • 



(2.12) 



with i,j running over 1,3 (or any other pair). To determine the coefficients of integrals 
we must sum over the two solutions corresponding to 7+ and 7_. It turns out that for the 
three-external-mass case, these solutions are related by taking t —>■ 1/t. In addition, when 
a corner of the triangle is massless, simpler forms of the solutions can be obtained. These 



issues will be discussed elsewhere 
massive case 



61|. A similar solution to eq. (I2.9P has been given in the 



65|. 



OPP [3JJ showed that after subtracting the known box contributions from the triple cut 
integrand, one is left with seven independent coefficients. One of these seven corresponds to 
the coefficient of the scalar triangle we seek, while the remaining six correspond to terms that 
integrate to zero. Evaluating the subtracted triple-cut integrand at seven selected kinematic 
points leads to a system of linear equations for these coefficients. As discussed in ref. [351 ] . 
however, numerical stability issues can arise from inverting this linear system of equations. 
The OPP approach of solving a system of equations is currentl y being implemented in 



36 



37 



38 



39j. In the 



numerical programs, with initial results reported in refs. Ill, 

n 

alternative approach of Forde [34|, the coefficient is instead extracted from the analytic 
behavior of the triple cut in the limit that the complex variable t becomes large. 
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We choose to use a hybrid of these approaches, subtracting box contributions from the 
triple cuts following OPP, but in a way that makes manifest the analytic properties in the 
complex variable t following Forde. The triple cut is, 

C 9 (t) = A^A^A^\ k=hity (2.13) 

Each of the box contributions to the triple cut (I2.13P contains a fourth Feynman propagator, 
1/7? (t) for some i ^ 1, 2, 3. Hence C${t) develops a pole in t whenever the inverse propagator 
vanishes, say 

as t-tf. (2.14) 

The pole locations t? and coefficients £f are determined from the form of lf(t), after inserting 
the triple-cut loop momentum parametrization (12. 9p . 

The residues at the poles also involve the coefficients d" of the i th box integral, evaluated 
on the two solutions a to the quadruple cuts. The d\ can be computed prior to the triangle 
calculation, and their contribution subtracted to form the difference, 

m ^ C3{t) ~fJrW(^)- (2 ' 15) 

Equation (I2.15P is slightly schematic, omitting a few subtleties that depend in part on how 
many of the triangle legs are massive. For example, in the three-mass case we should either 
sum over 7+ and 7_, or else make use of the t —>■ 1/t relation between the two triple-cut 



solutions to eliminate one of them 6l|. The main point is that proper subtraction of the 
box contributions removes all poles at finite values of t, so that T%(t) has poles only at t — 
and t = 00, as sketched in fig. HI Thus we can write, 

Ut)= £ c 3 tK (2.16) 

3=-P 

From eq. (12. 9p we see that the maximum power of t in eq. (12.161) . denoted by p, is equal 
to the maximum tensor rank encountered at the level of triangle integrals. In a generic 
renormalizable theory such as QCD, this value is p = 3. 

As explained in ref. {34], the desired coefficient of the triangle integral is given by Cq, 
which can be extracted by taking the limit t — > 00 and keeping only the t° contribution. 
This "Inf" operation can be applied to either C%(t) or the box-subtracted triple cut integrand 
T 3 (t), because the box contributions vanish as t — ► 00. In the language of OPP, the terms 
with j 7^ in eq. (12.161) correspond to terms that integrate to zero. 
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FIG. 4: After subtracting the box contributions to the triple cut, the t plane is free of all singularities 
except at t = and t = oo. We can extract the desired triangle coefficient by using a discrete 
Fourier projection, evaluating T 3 (i) at points indicated by the squares on the circle. 



We can also express the triangle coefficient using a contour integral around t = 0, 

1 r dt m „ . . 
co = ^-f T T 3 (t), (2.17) 

as depicted in fig. @J However, because of the special analytic form (12.161) of T 3 (t), it is much 
more efficient numerically to evaluate this contour integral by means of a discrete Fourier 
projection, 

c = ^—t T 3 (t e 2 ^ /(2p+1) ), (2.18) 

ZP + 1 }=-p 

where to is an arbitrary complex number. This projection removes the remaining coefficients 
Cfe, k 7^ 0. As it turns out, we do need the other coefficients in order to subtract out triangle 
contributions when evaluating bubble coefficients 6l|. We can obtain them from the same 
2p+ 1 evaluations of T 3 (t), by multiplying or dividing by factors of t before carrying out the 
Fourier sum, 



c k = j^— £ [toe 2 ^/(^+i)]- fe T 3 (t e 2 ^/(2p + i)) . (2 . 19) 



3=-P 

As we shall discuss in section [TV] the discrete Fourier projection provides excellent numerical 
stability. 



C. Bubble Coefficients 



Next consider the bubble coefficients. To parametrize the remaining degrees of free- 
dom left by the two-particle cuts shown in fig. EJ^b), we make use of a lightlike vector Ky 

13 



constructed from the external momentum and an arbitrary lightlike vector x M . The asso- 
ciated spinors are \Kf) and \x ± )- The normalization of x M = (x I 7 M \x~) /2 is determined 
by the constraint that K i • x — ^i/2, which ensures that 



(2.20) 



is lightlike. Note that this definition of K\ differs from the one f 1 2 . X j) in the triangle discus- 
sion, and is used exclusively for the two-particle cuts associated with the bubble coefficient. 
The cut conditions If = (i = 1,2) are solved by the momenta, 

1 



2 Kt + (y-\)(K^ 



x") + l (^ririx-) + ^^(xiTi^r), (2.21) 



with two free parameters y and t 3l|, [34 ] . 

In the two-particle cuts it is sometimes useful to restrict the cut loop momenta to be real. 
In this case, for Si = K\ > 0, the cut corresponds to a physical rescattering process. It is 
convenient to view the rescattering in the center-of-mass frame, in which K% = (y/Si, 0, 0, 0), 
the energies of the intermediate momenta li(y, t) are fixed to be y/Si/2, and the phase space 
can be parametrized alternatively by the polar and azimuthal angles 9 and for one of the 
two momenta, say l\. The relation between the two parametrizations is given by, 



• 2$ 

sin - 



t = - sin 9 e i(p . 
2 



(2.22) 



Then y is real and restricted to y G [0, 1], while t = Jy(l — y) e 1 ^ with G [0, 2tt). 

After subtracting box and triangle contributions from the two-particle cut under consid- 



eration 



31 



C 2 (y,t)=Af^A^' 



k=h(y,t) 



(2.23) 



l (i) ^(2) 

we are left with a tensorial expression B2(y,t) in terms of the loop momentum Zj, with 
maximal rank (p — 1). (In general, if the maximal rank of the triangle integrals is p, the 
maximal rank of bubble integrals is p — 1.) In terms of the parametrization (12.211) . B 2 (y, t) is 
a (p— l) th order polynomial expression in terms of the monomials (1/2 — y), t and y(l — y)/t. 



The bubble coefficient is then given by the integral 61], 

1 r 1 r fit 
b = — dyf ^B 2 {y,t). (2.24) 

ITXI JO J\t\ = y/y{l-y) t 

The factor of 1/t is a Jacobian for the change of variables f!2.22j) from (9,<fi) to (y,t). 
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As in the case of the triangle coefficients, the special analytic form of the subtracted 
two-particle cut -82(1/, t) allows the integral (12.241) to be evaluated efficiently using a discrete 
Fourier projection. Two observations are important here: the t integration projects B 2 (y,t) 
onto the terms independent of t, which are of maximal power (p — 1) in y; also, t he y 
integration amounts to replacing positive powers of y n by rational numbers l/(n + 1) 34]. 
Following similar logic as in the triangle case, we can extract the bubble coefficient with a 
double discrete Fourier projection on the subtracted two-particle cut, 

1 2(p-l)p-lp-l / 2nik/p\-n 

bo = j^-- E EE [ -T B 2 (y e 2 ™^, t e **7C*-D) , (2.25) 
(2p-l)p ~ jS^o n + 1 V J 

where yo and t are arbitrary complex constants. For the case p = 3, we use the fact that 
for f(y) — fo + fiy + /2I/ 2 , the desired combination f + fi/2 + / 2 /3 can be written as 
[/(0) + 3/(2/3)]/4. In this way it is possible to reduce the number of values of y required, 
from three in eq. (I2.25P to two: 

4 

20 



b o = hi E I** (°» *o e 2 ^ 75 ) + SB 2 (2/3, t e 



2-Kij/b 



(2.26) 

3=0 

One can also reduce the number of values of t sampled, from five down to three or four, 
using lower-order roots of unity (independently of how y is treated). In a similar fash- 
ion to eq. (I2.19p . higher-rank tensor bubble coefficients may be extracted by weighting the 
sum (I2.25P differently. (Such coefficients would feed into the calculation of tadpole coeffi- 
cients. They are not needed for the case of massless internal lines treated in this paper.) 

Due to the physical interpretation of the two-particle cut as a rescattering, with real 
intermediate momenta living on a sphere, an alternative projection formula from eqs. (I2.25P 
and (I2.26P may be found in terms of spherical harmonics Y^ m (8,<j)). To do so we change 
from the variables y and t to the spherical coordinates 9 and (f) via eq. (I2.22p . In these 
variables, the loop momentum (I2.2ip is linear in the spherical harmonics Y^ m with I = 1 and 
m — 0, ±1, because 

1 1 / 

2 -y = 2 cos ^ = V3*i,o(M) » 



t = Isintfe^-Wyli 



^ = jflWM- (2.27) 
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FIG. 5: Using Cauchy's theorem, rational terms in loop amplitudes can be reconstructed from 
residues at poles in the complex plane. The poles are of two types: physical and spurious. All 
pole locations are known a priori. Residues at physical poles are obtained from on-shell recursion. 
Residues at spurious poles are obtained from the cut parts. 



The two-particle cut with box and triangle contributions subtracted is then a superposition 
of spherical harmonics, 

B 2 (6A)= E h,mYi, m {9,<l>). (2.28) 
M<i<p-i 

The scalar bubble coefficient is just &o,o> U P to a normalization constant. Using eq. (I2.2ip . the 
higher spherical-harmonic coefficients h^ m can be related to the coefficients of the higher-rank 
tensor integrals. 



III. RATIONAL CONTRIBUTIONS 

We now turn to the question of computing the rational terms R n in the amplitude fl2.ll) . 



271 ]. modifying it 



Here we use the on-shell recursive approach for one-loop amplitudes |26 
to make it more amenable to numerical evaluation in an automated program. As is true for 
the cut parts, an important feature of on-shell recursion is that it displays a modest growth 
in computational resource requirements — compared to the rapid growth with a traditional 
Feynman-diagram approach — as the number of external particles increases. 

At one loop, as at tree level, on-shell recursion provides a systematic means of determining 
rational functions, using knowledge of their poles and residues. At loop level, however, a 
number of new issues must be addressed, including the appearance of branch cuts, spurious 
singularities, and the behavior of loop amplitudes under large complex deformations. In 
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some cases, "unreal poles" develop [251 ]. which are poles present with complex but not real 



momenta. The appearance of branch cuts does not present any difficulties because we use 



on-shell recursion only for the cut-free rational remainders R n . As noted in ref. [271 ]. we 
can sidestep the problems of unreal poles by choosing appropriate shifts within the class 
given below in eq. fl3.ll) . Finally, we may determine the behavior of amplitudes under large 
complex deformations by using auxiliary recursion relations. 

A. General Principles 

On-shell recursion relations may be derived by considering deformations of amplitudes 
characterized by a single complex parameter z, such that all external momenta are left on 
shell {24J]. In the massless case, it is particularly convenient to shift the momenta of two 
external legs, say j and /, 

e s -> q(z) = kf - !<ri7in , 

fcf ->TO = tf + fO-|7 M |l->- C 3 - 1 ) 

We denote the shift in eq. (13.11) as a [j, I) shift. This shift has the required property that 
the momentum conservation is left undisturbed, while shifted momenta are left on-shell, 
k](z) = kf(z) = 0. 

On-shell recursion relations follow from evaluating the contour integral, 

2m Jc z 

where the contour is taken around the circle at infinity, as depicted in fig. [5], and R n (z) is R n 
evaluated at the shifted momenta (I3.ip . If the rational terms under consideration vanish as 
z — > 00, the contour integral vanishes and Cauchy's theorem gives us a relationship between 
the desired rational contributions at z = 0, and a sum over residues of the poles of R n (z), 
located at z a , 

R n (0) = - £ Res^. (3.3) 

poles a 

On the other hand, if the amplitude does not vanish as z — > 00, there are additional contri- 



butions, which we can obtain from an auxiliary recursion relation [271 ]. 

Poles in the z-shifted one- loop rational terms, labeled by a in eq. (13. 3j) . may be separated 
into two classes as shown in fig. [5j physical and spurious. The physical poles are present in 
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FIG. 6: Diagrammatic contributions to on-shell recursion at one-loop for a [j, I) shift. The labels 
"T" and "L" refer respectively to (lower-point) tree amplitudes A tiee and rational parts of one- 
loop amplitudes R. The central blob in (c) is the rational part of a one-loop factorization function 



66] 



the full amplitude A n , and correspond to genuine, physical factorization poles (collinear or 
multiparticle). The spurious poles are not poles of A n ; they cancel between the cut parts 
C n and rational parts R n . They arise from the presence of tensor integrals in the underlying 
field-theory representation of the amplitude. Our method avoids the need to perform the 
reduction of such tensor integrals explicitly, because of the use of a basis of master integrals. 
The reduction happens implicitly, and leaves its trace in the presence of Gram determinant 
denominators. These denominators give rise to spurious singularities in individual terms. 
Separating the different contributions, we may write, 

R n (z) = R°{z) + R s n {z) + < argc *(*) , (3-4) 

where R^ contains all contributions from physical poles, R^ the contributions from spurious 
poles, and R l £ rge z the possible contributions from large deformation parameter z, if R n {z) 
does not vanish there. More explicitly, from elementary complex variable theory, under the 
shift (13. ip the rational terms can be expressed as a sum over pole terms and possibly a 
polynomial in z, 

R D niz) = E — > = E { + — ) ' 

a Z~Z a y\{Z-ZpY Z-Zp) 

& max 

RT eez (z)= E AX, (3.5) 

(7=0 

where the coefficients A a , Bp, Cp, D a are functions of the external momenta. The poles in z 
in eq. (13. 5p are shown in fig. The physical poles labeled by a are generically single poles. 
(Some shift choices may lead to double poles j^J; we can generally avoid such shifts ^\.) 
In general, in a renormalizable gauge theory, the spurious poles, labeled by /3, may be either 
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single or double poles [6l|]. If R n (z) vanishes for large z, the D a are all zero. If not, then 
Do gives a contribution to the physical rational terms, R n {0). 

The contributions of the physical poles may be obtained efficiently using the on-shell 
recursive terms represented by the diagrams in fig. [6J The tree "vertices" labeled by "T" 
denote tree-level on-shell amplitudes while the loop vertices "L"are the rational parts 

of on-shell (lower-point) one-loop amplitudes R m , m < n, as defined in eq. (I2.3p . The 
contribution in fig. |6(c) involves the rational part of the additional factorization function 
T 0]. It only appears in multi-particle channels, and only if the tree amplitude contains 
a pole in that channel. Each diagram is associated with a physical pole in the z plane, 
illustrated in fig. 0, whose location is given by, 



(3.6) 



where K- 



K 



r...s\^rsj 



(j \$ r ... s \l ) 

= k r + k r+1 + • - • + k 8 . This pole arises from the vanishing of shifted propagators, 
0. Generically the sum over a is replaced by a double sum over r, s, labeling 



the recursive diagrams, where legs labeled j and / always appear on opposite sides of the 
propagator in fig. [61 The computation of the recursive diagrams has been described in 



refs. 



26 



33 



44|, to which we refer the reader for further details. 



What about the contributions of the spurious poles? One approach is to find a "cut 



completion" 



26 



271 ]. which is designed by adding appropriate rational terms to C n in order 



to cancel entirely the spurious poles in z within the redefined cut terms C n . Because the 
complete amplitude is free of the spurious poles, this procedure ensures that the redefined 
rational terms R n are free of them. The cut completion makes it unnecessary to compute 
residues of spurious poles (although additional "overlap" diagrams are introduced). It is 
very helpful for deriving compact analytic expressions for the amplitudes. This approach 



has led to the computation of the rational terms 
with an arbitrary number of external legs 



26 



27 



br a variety of one-loop MHV amplitudes 
as well as for six-point amplitudes. In 



general, it should be possible to construct a set of cut completions using integral functions 
of the type given in ref. 54( to absorb spurious singularities. 

For the purposes of a numerical program, however, it is simpler to extract the spurious 
residues from the known cut parts. These residues are guaranteed to be the negatives of the 
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spurious-pole residues in the rational part. That is, the spurious contributions are, 

R s n (0) = - V Res^^= V Res^^, (3.7) 

^ z=Zg Z z=z Z 

spur, poles (3 " spur, poles f3 

where C n (z) is the shifted cut part appearing in eq. (12.1 1) . The spurious poles (3 correspond 
to the vanishing of shifted Gram determinants, A m (z) = for m = 2,3,4, associated with 
bubble, triangle and box integrals. (In the case of massless internal propagators, the bubble 
Gram determinant does not generate any spurious poles.) 

A simple example of a spurious singularity in the cut part (12.21) is from a bubble term of 
the form, 

hil= {K 2 h l K 2 2? H- K l) + ---> (3-8) 

where hi is smooth as K\ — > Iff, and Ki + K 2 + k 3 = for some massless momentum 
k%. The denominator factor (Kf — ff|) is the square root of the Gram determinant for a 
triangle integral with two massive legs, K\ and K 2 , and one massless leg, k 3 . Under the 
[j, I) shift, there will be a value of z, zp, for which the shifted denominator vanishes linearly, 
Kf(z) — K^{z) ~ z — Z/3 (unless j and I both belong to the same massive momentum cluster, 
K\ or K 2 , in which case the Gram determinant is unshifted). From eq. (13. 7p we see that we 
only need the rational pieces of the spurious-pole residues of the cut part, because (0) 
is rational. From eq. (13.81) . we see that there can only be a rational piece if we have to 
series expand the logarithm to compute the residue. Hence the spurious pole in the bubble 
coefficient bi must be of at least second order in (Kf — K%). At order e°, box and triangle 
integrals contain dilogarithms and squared logarithms, which must be expanded to second 
order to obtain a rational piece. Thus the spurious poles of box and triangle coefficients 
must be at least of third order for rational terms to be generated. 

To extract a residue from C n (z)/z, we evaluate the integral coefficients di,Ci,bi numeri- 
cally for complex, shifted momenta in the vicinity of the spurious pole, using our implemen- 
tation of the results of section [Til We also need to evaluate the loop integrals. First, however, 
we perform an analytic series expansion of the integrals around the vanishing Gram deter- 
minants. For example, the three-mass triangle integral, J| m (si, s 2 , s 3 ), close to the surface 
of its vanishing Gram determinant, 

A 3 = si + s\ + s\ - 2 Sl s 2 - 2sis 3 - 2s 2 s 3 -> 0, (3.9) 
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behaves as, 



1 3 

--^ln(- Si 

z i=l 



+ 



1 A a 



Si Si— 1 

s i+lSj-l 

Si + S 2 + S 3 / Aa 



1 A 3 1 
1 — + — 

6s i+ iSi_i 30 



A, 



+ 



(3.10) 



6sis 2 s 3 120 \sis 2 s 3 , 

where the index i on the shifted invariant, = Si(z), is defined mod 3. In this expression 
the logarithms are to be expanded according to, 

1 (s - spf 



ln(-s) 



s-sp 



+ 



(3.11) 



Sp I S*p 

where s = s(z), and Sp = s(zp) is the value of the shifted invariant at the location z a of the 
spurious pole. The leading order of eq. f 1 3 . 1 f) matches the expansion found in ref. j54|. In 
the integral expansions we need keep only rational terms, including terms that can become 
rational after further series expansion around a generic point, such as eq. (13.111) . Thus we 
may avoid computing any logarithms or polylogarithms at complex momentum values. The 
expression obtained by repla cing C n (z) according to these rules, in the vicinity of zp, will 
be denoted by E@(z). In ref. [6l| we present the complete set of integral expansions needed 
in the calculations, as well as a convenient method for generating them from a dimension- 
shifting formula 47]. 



B. Discrete Fourier Sum for Spurious Residues 

Similarly to the case of triangle and bubble coefficients, we extract each required spurious- 
pole residue from the cut parts by using a discrete Fourier sum. We evaluate E@(z) at m 
points equally spaced around a circle of radius 5p in the z plane, centered on the pole 
location zp, as depicted in fig. i.e., z = zp + 5pe 2m ^ m , for j — 1, 2, . . . , m. In contrast to 
the t-plane analysis used earlier to obtain triangle coefficients, however, we do not know the 
residues at other poles a priori, so we cannot subtract them easily. (Indeed, the function 
E@(z) we are analyzing is only rational in the vicinity of zp, due to our use of the rational 
parts of the integral expansions around this point.) Here the discrete Fourier sum is an 
approximation to the contour integral, whereas in the previous section it was exact. We can 
make the approximation arbitrarily accurate in principle, by choosing 5p to be arbitrarily 
small. With finite precision, however, numerical round-off error forces us to work at finite 
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\_z_ 




FIG. 7: We obtain the residue at the spurious pole located at z = zp in the complex z plane by 
a discrete Fourier sum, evaluating E@(z) on the (blue) squares on the circle of radius 5p centered 
on zp. In this figure four points are shown, although in practice we use ten points. The locations 
of other poles are represented by (red) dots. We ensure that 5p is sufficiently small so that other 
poles give a negligible contribution to the residue. 



5 p. When extracting the residue of a spurious pole we must also ensure that there are no 
other poles inside or near the circle. To obtain the contributions of the spurious poles to 
Rn(0) in eq. (13. 7p we evaluate, 

rrippl zp + 6pe 2m ^ m 

The sum over (3 runs over the location of all spurious Gram determinant poles that contribute 
to rational terms. Equivalently, we can extract the coefficients Bp and Cp in eq. ( 13. 51) via, 



-I m 2 

Bp ~ Y\5 p e 2 ^' m ] E^(zp + Spe 2m ^ m ), 

m 3=1 
i m 

Cp ~ J2 5 P e27Tij/m E n( z P + 6pe 2%ij/m ) . (3. 13) 

m 3=1 

For the results presented in the next section we choose m = 10 points in the discrete sum. 
In general, an increase in m increases the precision, but at the cost of computation time. 

We choose dp to be much smaller than the distance to nearby poles, but not so small 
as to lose numerical precision. Typically at "standard" double precision we use a value of 
5p = 10~ 2 . If the contributions from the nearby poles are unusually large, then we find a 
large variation in the absolute value of each term in the sum. If this happens we reduce 5p 
until either the variation is acceptable, or we cross a minimum value of Sp, beyond which the 
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point becomes unstable because of round-off error. We deal with such points as described 
below. 



C. Numerical Stability 

In addition to the value of 5p becoming too small, other cancellations can also sometimes 
cause a loss of precision, giving rise to a potentially unstable kinematic point. In order to 
identify such phase-space points more generally, we apply consistency checks independently 
to the cut and rational parts of the amplitude. For the cut part we test how well the known, 
non- logarithmic 1/e singularities are reproduced. Because the only source of such 1/e poles 
are the bubble integrals, for the n-gluon amplitudes, for example, we have [a, 1101]. 

<- loop |iA,„on^o g = Ip* = " [7 (y - |^)] AT, (3.14) 

where rif is the number of quark flavors and the sum on k runs over all bubble integrals. As a 
practical matter it is sufficient to check that the divergent term divided by the tree amplitude 
is real. (For helicity configurations with vanishing tree amplitudes the cut contributions 
vanish, so no check is required.) Because bubble coefficients are computed from expressions 
where triangle and box contributions have been subtracted, any instabilities in the latter 
are also detected with this 1/e consistency check. 

In general this test is not sufficient for finding all the unstable points of the full amplitude, 
because some of the instability comes from computing the spurious residues for rational 
terms. A related test, which suffices to find all remaining instabilities, comes from the 
requirement that each spurious singularity must cancel in the sum over bubble coefficients. 
This cancellation can be understood by applying the [j, I) shift to eq. (13.141) . and making use 
of the fact that A^ ee has no spurious poles. For each spurious-pole residue that contributes 
to the rational part, we therefore check that the sum of discrete Fourier sums over all bubble 
coefficients, 

m 

E E h^ lm h{z P + V 2 " j/m ) > (3-15) 

k 3=1 

vanishes to within a specified tolerance. 

If a phase-space point fails the above stability conditions we recalculate the point in a 
manner that improves its stability. Various strategies have been proposed in the literature 
to handle unstable points. One approach is to modify the standard integral basis (12. 2p so as 
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to absorb the Gram determinant singularities into well-defined functions 



54 



67|- 



26J. Other approaches are to interpolate 



This approach is related to using a cut completion 

across the singular region or to series expand the integrals in the singular region {5 
third approach is to simply redo unstable points at higher precision, e.g. as in ref. 36] 
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We have found the high-precision approach to be effective for eliminating the remaining 
instabilities in our program. It is robust and simple to implement; a detailed analysis 
of the instabilities is not needed, and we can use the standard basis of integrals with no 
interpolations or expansions of the integrals around unstable points. Our implementation 
of on-shell methods already has only a small fraction of unstable phase-space points; hence 
the overhead of recomputing them at higher precision is relatively small. We use the QD 
package [5^, switching to "double-double" precision, that is approximately 32 decimal digits. 
If the stability test were to fail at this level of precision, we switch to "quadruple-double" 
precision, corresponding to approximately 64 digits of precision; for all amplitudes calculated 
here, this happens rarely, if ever. To compute the integrals to higher precision, we implement 
the polylogarithms which enter the integrals using a series expansion to a sufficiently high 
order. If the 1/e test (13. 14ft fails then we recompute the entire cut part at higher precision, 
but if the spurious-pole test (13.151) fails we only recompute those pieces containing unstable 
Gram determinant singularities. 

Further details, as well as all integral expansions used to extract the spurious residues 
from the cut part, will be given elsewhere 61|. 



IV. RESULTS 



We now discuss the numerical stability of our implementation. Our stability tests use 
sets of 100,000 points for 2 — > (n — 2) gluon scattering, generated with a flat phase-space 
distribution using the RAMB0 68] algorithm. We impose kinematic cuts on the outgoing 



gluons, following ref. 



35j: 



E T > 0.0lVs ; 



V <3 



An > 0.4 



(4.1) 



where E T is the gluon transverse energy, 77 is the pseudorapidity, and A^ = y A,^ + is 
the separation cut between pairs of gluons. The center-of-mass energy y/s is chosen to be 2 
TeV and the scale parameter fi (arising from divergent loop integrals) is set to 1 TeV. 
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FIG. 8: The distribution of the logarithm of the relative error for 100,000 phase-space points in the 
1/e 2 , 1/e and finite (e ) components of the six-point MHV amplitude Aq(1~ , 2~ , 3 + , 4 + , 5 + , 6 + ). 
The solid (black) curve shows the distribution run entirely with ordinary double precision, and the 
dashed (red) curve shows it when contributions identified as unstable — following the discussion 
of section jTnl — are evaluated using higher precision. The target values use analytic results from 



refs. 



13, 



42] 



We computed one-loop six-, seven- and eight-gluon amplitudes for rif = with BlackHat 
at each phase-space point, and compared the output against a target expression, obtained 
either from known analytic results, or from BlackHat itself using quadruple-double preci- 
sion (~64 digits). As an additional test, we also used ordinary double precision to compare 



to the numerical results of refs. 
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37] at the quoted phase-space points. We find agree- 



ment for the five- and six-gluon amplitudes for all helicity configurations, to within their 
quoted accuracy, after accounting for external phase conventions and the incoming-particle 
convention implicitly used in ref. 37]]. We also find agreement with the numerical results 



of ref. 



27] at the quoted phase-space points for the six-, seven- and eight-point maximally 



helicity violating (MHV) amplitudes presented here. 

The histograms in figs. IBliTOl show the results of our study of numerical precision. For 
these plots, the horizontal axis is the logarithmic relative error, 

I yjnum ^target I \ 



log 



10 



Kin 



target i 



(4.2) 



for each of the 1/e 2 , 1/e and e° components of the numerical amplitude A™ m obtained from 
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BlackHat. The vertical axis in these plots shows the number of phase-space points in a bin 
that agree with the target to a specified relative precision. We use a logarithmic vertical scale 
to visually enhance the tail of the distribution, so as to illustrate the numerical stability. 
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For the MHV amplitudes in figs. [S] and [HI we used analytic expressions from refs. 



13, 



421 . |44j as the target expressions A^* get . 



analytic expressions are available 
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n 



29 



or the next-to-MHV (NMHV) amplitudes, 



42j | . although for fig. [TUl we generated the 



target with BlackHat, using quadruple-double precision. This is more than sufficient to 
ensure numerical stability in target expressions for the purposes of the comparison. We note 
that the ability to switch easily to higher precision is quite helpful in assessing numerical 
stability in any new calculation. 

First consider the MHV six-point amplitude A$(l~, 2~, 3 + , 4 + , 5 + , 6 + ). Fig. [H] illustrates 
the numerical stability of BlackHat for this amplitude, with and without the use of higher 
precision on the points identified as unstable. The plots show the distribution of relative 
errors for the 1/e 2 , 1/e and e° components over 100,000 phase-space points. The 1/e 2 
distribution has extremely small errors, peaking at a relative error of nearly 10 -15 , while 
the right-side tail falls rapidly. For the 1/e and finite e° components the peaks shift to the 
right, to a relative precision of around 10~ 14 and 10~ n , and fall less steeply. This feature is 
not surprising, because of the larger number of computational steps needed for these parts 
of the amplitudes: for 1/e 2 terms, only box coefficients contribute (for this helicity pattern 
triangle integrals do not appear); for the 1/e contribution, bubble coefficients contribute too; 
for the finite part, rational terms contribute as well. As one proceeds from box to triangle, 
bubble, and then to rational terms, each step relies on previous steps, and so numerical 
errors accumulate. 

In each plot in fig. [8] the solid (black) curve corresponds to the exclusive use of ordinary 
double precision (16 decimal digits), showing good stability for the raw algorithm for all 
three components. The dashed (red) curve shows the effect of turning on higher precision 
for contributions identified as unstable, using the criteria discussed in section IIHI This 
completely suppresses the already-small tail above a relative error of about 10~ 5 . The 
points populating the right-hand tail in the ordinary double precision calculation, displayed 
in the solid (black) curve, then move to the left in the dashed (red) curve, giving rise to 
a secondary peak around a relative error of machine precision, or 10 -16 . (The comparison 
with the target is performed in ordinary double-precision, even though higher precision is 
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FIG. 9: The distribution of the logarithm of the relative error over 100,000 phase-space 
points for the MHV amplitudes A & (1~ , 2~ , 3+, 4+, 5+, 6+), A 7 (l-, 2~, 3+, 4+ 5 5+, 6+, 7+) and 
As(l~ , 2~, 3 + , 4 + , 5 + , 6 + , 7 + , 8 + ). The dashed (black) curve in each histogram gives the relative 
error for the 1/e 2 part, the solid (red) curve gives the 1/e singularity, and the shaded (blue) distri- 
bution gives the finite e° component of the corresponding helicity amplitude. The target expression 



is computed from an analytic formula 



13, 



26 



42 



44J. 



used in intermediate steps.) This twin-peak feature is visible in the 1/e and e° components. 
It is due our recalculation of the entire cut part, at higher precision, whenever a phase-space 
point fails the 1/e consistency check (13.141) . When the spurious-pole stability test (13 . 15[) 
fails, the point generally falls to the right of the secondary peak, because we only recalculate 
those pieces that contain the unstable spurious singularity. 

Another important feature that can be observed in fig. [8] is that the "effective cutoff" 
is sharp: for the e° terms almost no points below 10~ 5 are identified as unstable. In a 
practical calculation, given Monte-Carlo integration errors and other uncertainties, a cutoff 
in the relative error of 10~ 5 is overly stringent. It does, however, illustrate the control over 
instabilities achieved in BlackHat, which becomes more important for more complicated 
processes. It is interesting to note that modest additional computation time is required to 
achieve a cutoff of 10~ 5 , compared to, say, 1CT 2 . 

Next consider the behavior as the number of external gluons increases. In fig. [9] we 
show relative error distributions for the set of MHV amplitudes Aq(1~, 2~, 3 + , 4 + , 5 + , 6 + ), 
A 7 (l-, 2", 3+ 4+ 5+ 6+, 7+) and A 8 (l~, 2" 3+ 4+ 5+ 6+ 7 + , 8+). For each of these ampli- 
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FIG. 10: The distribution of the logarithm of the relative error for the six-point NMHV amplitudes 
A;(l-,2-,3-,4 + ,5 + ,6 + ), ^ 6 (l-,2-,3 + ,4-,5+,6+) and A 6 (l~, 2+, 3", 4+, 5~, 6+). The dashed 
(black) curve in each histogram gives the relative error for the 1/e 2 part, the solid (red) curve gives 
the 1/e singularity, and the shaded (blue) curve gives the finite e° component of the corresponding 
amplitude. The target expression is a quadruple-double-precision BlackHat evaluation. 

tudes the dashed (black) curve shows the relative error in the coefficient of the 1/e 2 singular- 
ity. Similarly, the relative errors in the 1/e and e° contributions are given by the solid (red) 
curve and shaded (blue) distribution. The relative precision of the 1/e 2 singularities is better 
than 10~ n for these six-, seven- and eight-point amplitudes. The computational-stability 
scaling properties in going from six- to seven- and then eight-point amplitudes in fig. [9] are 
also rather striking. There is little change in the shape of the curves as we increase the 
number of legs. 

Even more striking is the modest increase in computation time. As mentioned earlier, the 
tree-like nature of on-shell methods leads us to expect only mild scaling for a given helicity 
pattern, in stark contrast with the rapid increase in required computational resources for 
ordinary Feynman diagrams. These expectations are borne out by the values for the average 
computation time shown in Table [H The table shows the average time on a 2.33 GHz Xeon 
processor for computing a color-ordered amplitude of a given helicity configuration at a 
single phase-space point. The first three rows show the timing for the six-, seven- and eight- 
point MHV amplitudes corresponding to fig. [9j Even for the eight-point case we obtain an 
average evaluation time of less than 50 ms, including running the phase-space points marked 
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TABLE I: The average time needed to evaluate one point in phase space for various helicity 
configurations. The time is in milliseconds on a 2.33 GHz Xeon processor. The second column 
gives the average evaluation time for the cut part, including the recomputation at higher precision 
of points identified as unstable. The third column gives the time for the full amplitude, including 
rational terms, using only ordinary double precision. The fourth column gives the average time 
using ordinary double precision on stable points and higher precision on contributions marked as 
unstable either by the 1/e consistency test (|3.14p or the spurious-pole test ()3. 15j) . 



helicity 


cut part 


full amplitude 
double prec. only 


full amplitude 
with multi-prec. 


++++ 


2.4 ms 


7 ms 


11 ms 




4.2 ms 


11 ms 


23 ms 


++++++ 


6.1 ms 


29 ms 


43 ms 


-+-+++ 


3.1 ms 


18 ms 


32 ms 


-++-++ 


3.3 ms 


61 ms 


96 ms 


+++ 


4.4 ms 


12 ms 


22 ms 


+-++ 


5.9 ms 


47 ms 


64 ms 


-+-+-+ 


7.0 ms 


72 ms 


114 ms 



as unstable at higher precision. It is also interesting to note the relatively modest increase 
in computation time due to turning on higher precision for unstable points, even in this 
initial implementation. (The time in the third column includes the evaluation of bubble 
coefficients used in the spurious-pole test ( 13.151) .) 

Finally, consider the six-gluon NMHV amplitudes. Figure [10] illustrates the numer- 
ical stability properties of the complete set of independent six-gluon NMHV ampli- 
tudes not related by symmetries, A 6 (l~ , 2~, 3~, 4 + , 5 + , 6 + ), A e (l~ , 2~, 3 + , 4~, 5 + , 6 + ) and 
Aq(1~ , 2 + , 3~, 4 + , 5~, 6 + ), compared against a quadruple-precision target computed with 
BlackHat. For each one of these amplitudes, the contributions to the 1/e 2 , 1/e and fi- 
nite e° terms are shown in a similar format as the MHV case. These NMHV curves are all 
shifted to the right compared to the MHV cases in in fig. [9j This property is not surprising; 
it is due to the more complicated nature of the NMHV amplitudes. In particular, the am- 
plitudes contain higher powers of the box Gram determinants in denominators of the box 
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coefficients, which then feed into triangle, bubble and rational contributions. As in the MHV 
cases, when one goes from 1/e 2 to 1/e to e°, the curves shift to the right again, reflecting the 
more complicated calculations. Nevertheless, they all exhibit excellent numerical stability, 
with the distributions of relative errors for the finite pieces peaking at 10~ 8 or better. We 
identify points as unstable, and automatically recompute such points at higher precision, 
using the same criteria as for the MHV amplitudes. In the NMHV case, the fall-off is not 
as sharp as in the MHV case. Nevertheless, the accuracy obtained is more than sufficient 
for use in an NLO program. 

The average evaluation time in the current version, for all independent six-gluon helicity 
configurations needed at NLO, including the NMHV ones, is given in Table [B One can see 
that alternating-helicity configurations do take longer to compute. However, in all cases the 
cut parts are evaluated in under 8 ms and the full amplitudes in under 120 ms. Although 
we have not run systematic tests of NMHV amplitudes beyond six points, initial studies 
at seven points indicate that the scaling behavior of the NMHV amplitudes is not quite as 
good as for the MHV case, but still very good. 



V. CONCLUSIONS 



In this paper we presented the first results from BlackHat, an automated implementation 
of on-shell methods, focusing on the key practical issues of numerical stability and computa- 
tional time. We illustrated the numerical stability by computing a variety of complete six-, 
seven- and eight-gluon helicity amplitudes and comparing the results against previously- 
obtained analytic results or against higher precision calculations. In this initial version we 
achieved reasonable speed, an average computation time of 114 ms per phase-space point 
for the most complicated of the six-gluon helicity amplitudes, and substantially better for 
the simpler helicities. We expect this speed and stability to be sufficient for carrying out 
phenomenological studies of backgrounds at the LHC, even as we expect further improve- 
ments with continuing optimization of the code. After the code is stable and tested for a 
wide variety of processes, we plan to make it publicly available. 

BlackHat uses the unitarity method with four- dimensional loop momenta 13, |42j. This 
method allows the use of compact tree-level helicity amplitudes as the basic building blocks. 
We compute the box coefficients using quadruple cuts [21| . For box integrals with massless 
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internal propagators and at least one massless corner, we presented a simple solution to the 
cut conditions. The solution makes manifest the presence of square roots, rather than full 
powers, of a spurious (Gram determinant) singularity for each power of the loop momentum 
in the numerator. We evaluated the triangle- and bubble-integral coefficients using Forde's 
approach 34( to expose their complex- analytic structure. Another important ingredient in 
our procedure is the OPP [31] subtraction of boxes from triple cuts when computing triangle 
coefficients, and of boxes and triangles from ordinary (double) cuts when computing bubble 
coefficients. Viewed in terms of Forde's complex-valued parametrization approach, the OPP 
subtraction cleans the complex plane of poles, using previously-computed coefficients. We 
then introduced a discrete Fourier projection, as an efficient and numerically stable method 
for extracting the desired coefficients. In the bubble case, this procedure can be recast in 
terms of spherical harmonics. 

We computed the purely rational terms using loop-level on-shell recursion, modifying the 
treatment of spurious singularities compared to refs. [26|, [27|. We used a discrete Fourier 
sum to compute the spurious-pole residues from the cut parts. These contributions are 
then subtracted from the recursively-computed rational terms in order to cancel spurious 
singularities implicit in the latter, and thereby make the full amplitude free of spurious 
singularities as required. 

The computation of most points in phase space proceeds using ordinary double-precision 
arithmetic to an accuracy of 10~ 5 or less. This is far better than the Monte-Carlo inte- 
gration errors that will inevitably arise in any use of amplitudes in an NLO parton-level 
or parton-shower code (not to mention parton distribution, scale, shower and hadroniza- 
tion uncertainties). Nonetheless, the computation of the amplitude at a small percentage 
of phase-space points does manifest a loss of precision, resulting in an instability and larger 
error. In order to identify such unstable points as may arise, we impose the requirements 
that all spurious singularities cancel amongst bubble coefficients, and that the coefficients of 
the 1/e singularity (corresponding to e-singular terms in bubble integrals) be correct. When- 
ever the calculation at a given phase-space point fails these criteria we simply recalculate the 
point at higher precision. There are other possible means for dealing with Gram-determinant 
singularities Q, 55, 56 ] , but we prefer this approach because of its simplicity j^f]]. In 



practice, it has a relatively modest impact on the overall speed of the program. In the most 
complicated of the six-gluon helicity amplitudes, higher-precision evaluation causes the time 
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to increase modestly, from 72 ms to 114 ms. We expect to see further improvements with 
additional refinements. 

It is important to validate a numerical method against known analytic results. For 
this purpose, we made use of MHV configurations, which contain two gluons of helicity 
opposite to that of the others. In particular, we considered the case where the two opposite 



helicities are nearest neighbors in the color order. In ear 



ier work, these amplitudes were 



computed for an arbitrary number of external gluons [26|, |44j, using on-shell methods. We 
used these results to confirm that BlackHat returns the correct values through eight gluons. 
We also verified numerical stability for non-MHV amplitudes by comparing results for all six- 
gluon amplitudes against a reference computation done entirely using quadruple-precision 
arithmetic. 

We defer discussion of amplitudes with external fermions, or with massive quarks and 
vector bosons, to the future. (Some work directly relevant to the question of adding massive 



particles may be found in refs. [32j, |65|, |69(.) We will also present further details, including 
the integral expansions we use around spurious singularities, in a future publication 6ll |. 

The excellent numerical stability and timing performance of BlackHat is due to a variety 
of ideas described in this paper. Because the unitarity method uses gauge-invariant tree 
amplitudes as the basic input into the calculation, we avoid the large gauge cancellations 
inherent in Feynman- diagram calculations. In addition we made use of very compact four- 
dimensional tree-level helicity amplitudes as the basic input to the calculations. All steps 
in our computation of the rational terms, as well as the integral coefficients, are carried out 
in four dimensions. Our simple quadruple-cut solution (12. 7\\ also helps maintain numeri- 
cal stability in the box contributions. Our parametrization choices for triple and double 
cuts, and the OPP subtraction of previously-computed coefficients are additional important 
ingredients. Finally, our use of discrete Fourier projections helps considerably. 

The resulting C++ code BlackHat is efficient and numerically stable, as we have illus- 
trated with the computation of various one-loop gluon amplitudes and their comparison to 
known analytic expressions. Based on the results presented here, we expect BlackHat to 
make possible the computation of a wide variety of new one-loop amplitudes for collider 
physics that have been inaccessible with traditional methods. We hope that BlackHat, in 
conjunction with automated programs [9( for combining the real and virtual contributions 
at NLO, will soon enable the computation of phenomeno logically important cross-sections 
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at the LHC. 
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